Nonuniformity of P-values Can Occur Early in Diverging Dimensions

Evaluating the joint significance of covariates is of fundamental importance in a wide range of applications. To this end, p-values are frequently employed and produced by algorithms that are powered by classical large-sample asymptotic theory. It is well known that the conventional p-values in Gaussian linear model are valid even when the dimensionality is a non-vanishing fraction of the sample size, but can break down when the design matrix becomes singular in higher dimensions or when the error distribution deviates from Gaussianity. A natural question is when the conventional p-values in generalized linear models become invalid in diverging dimensions. We establish that such a breakdown can occur early in nonlinear models. Our theoretical characterizations are confirmed by simulation studies.


Introduction
In many applications it is often desirable to evaluate the significance of covariates in a predictive model for some response of interest. Identifying a set of significant covariates can facilitate domain experts to further probe their causal relationships with the response. Ruling out insignificant covariates can also help reduce the fraction of false discoveries and narrow down the scope of follow-up experimental studies by scientists. These tasks certainly require an accurate measure of feature significance in finite samples. The tool of p-values has provided a powerful framework for such investigations.
As p-values are routinely produced by algorithms, practitioners should perhaps be aware that those p-values are usually based on classical large-sample asymptotic theory. For ex-ample, marginal p-values have been employed frequently in large-scale applications when the number of covariates p greatly exceeds the number of observations n. Those p-values are based on marginal regression models linking each individual covariate to the response separately. In these marginal regression models, the ratio of sample size to model dimensionality is equal to n, which results in justified p-values as sample size increases. Yet due to the correlations among the covariates, we often would like to investigate the joint significance of a covariate in a regression model conditional on all other covariates, which is the main focus of this paper. A natural question is whether conventional joint p-values continue to be valid in the regime of diverging dimensionality p.
It is well known that fitting the linear regression model with p > n using the ordinary least squares can lead to perfect fit giving rise to zero residual vector, which renders the p-values undefined. When p ≤ n and the design matrix is nonsingular, the p-values in the linear regression model are well defined and valid thanks to the exact normality of the least-squares estimator when the random error is Gaussian and the design matrix is deterministic. When the error is non-Gaussian, Huber (1973) showed that the least-squares estimator can still be asymptotically normal under the assumption of p = o(n), but is generally no longer normal when p = o(n) fails to hold, making the conventional p-values inaccurate in higher dimensions. For the asymptotic properties of M -estimators for robust regression, see, for example, Huber (1973); Portnoy (1984Portnoy ( , 1985 for the case of diverging dimensionality p = o(n) and ;  for the scenario when the dimensionality p grows proportionally to sample size n.
We have seen that the conventional p-values for the least-squares estimator in linear regression model can start behaving wildly and become invalid when the dimensionality p is of the same order as sample size n and the error distribution deviates from Gaussianity. A natural question is whether similar phenomenon holds for the conventional p-values for the maximum likelihood estimator (MLE) in the setting of diverging-dimensional nonlinear models. More specifically, we aim to answer the question of whether p ∼ n is still the breakdown point of the conventional p-values when we move away from the regime of linear regression model, where ∼ stands for asymptotic order. To simplify the technical presentation, in this paper we adopt the generalized linear model (GLM) as a specific family of nonlinear models (McCullagh and Nelder, 1989). The GLM with a canonical link assumes that the conditional distribution of y given X belongs to the canonical exponential family, having the following density function with respect to some fixed measure where X = (x 1 , · · · , x p ) is an n × p design matrix with x j = (x 1j , · · · , x nj ) T , j = 1, · · · , p, y = (y 1 , · · · , y n ) T is an n-dimensional response vector, β = (β 1 , · · · , β p ) T is a p-dimensional regression coefficient vector, {f 0 (y; θ) : θ ∈ R} is a family of distributions in the regular exponential family with dispersion parameter φ ∈ (0, ∞), and θ = (θ 1 , · · · , θ n ) T = Xβ. As is common in GLM, the function b(θ) in (1) is implicitly assumed to be twice continuously differentiable with b (θ) always positive. Popularly used GLMs include the linear regression model, logistic regression model, and Poisson regression model for continuous, binary, and count data of responses, respectively.
The key innovation of our paper is the formal justification that the conventional pvalues in nonlinear models of GLMs can become invalid in diverging dimensions and such a breakdown can occur much earlier than in linear models, which spells out a fundamental difference between linear models and nonlinear models. To begin investigating p-values in diverging-dimensional GLMs, let us gain some insights into this problem by looking at the specific case of logistic regression. Recently, Candès (2016) established an interesting phase transition phenomenon of perfect hyperplane separation for high-dimensional classification with an elegant probabilistic argument. Suppose we are given a random design matrix X ∼ N (0, I n ⊗ I p ) and arbitrary binary y i 's that are not all the same. The phase transition of perfect hyperplane separation happens at the point p/n = 1/2. With such a separating hyperplane, there exist some β * ∈ R p and t ∈ R such that x T i β * > t for all cases y i = 1 and x T i β * < t for all controls y i = 0. Let us fit a logistic regression model with an intercept. It is easy to show that multiplying the vector (−t, (β * ) T ) T by a divergence sequence of positive numbers c, we can obtain a sequence of logistic regression fits with the fitted response vector approaching y = (y 1 , · · · , y n ) T as c → ∞. As a consequence, the MLE algorithm can return a pretty wild estimate that is close to infinity in topology when the algorithm is set to stop. Clearly, in such a case the p-value of the MLE is no longer justified and meaningful. The results in Candès (2016) have two important implications. First, such results reveal that unlike in linear models, p-values in nonlinear models can break down and behave wildly when p/n is of order 1/2; see ;  and discussions below. Second, these results motivate us to characterize the breakdown point of p-values in nonlinear GLMs with p ∼ n α 0 in the regime of α 0 ∈ [0, 1/2). In fact, our results show that the breakdown point can be even much earlier than n/2.
It is worth mentioning that our work is different in goals from the limited but growing literature on p-values for high-dimensional nonlinear models, and makes novel contributions to such a problem. The key distinction is that existing work has focused primarily on identifying the scenarios in which conventional p-values or their modifications continue to be valid with some sparsity assumption limiting the growth of intrinsic dimensions. For example, Fan and Peng (2004) established the oracle property including the asymptotic normality for nonconcave penalized likelihood estimators in the scenario of p = o(n 1/5 ), while Fan and Lv (2011) extended their results to the GLM setting of non-polynomial (NP) dimensionality. In the latter work, the p-values were proved to be valid under the assumption that the intrinsic dimensionality s = o(n 1/3 ). More recent work on high-dimensional inference in nonlinear model settings includes van de Geer et al. (2014); Athey et al. (2016) under sparsity assumptions. In addition, two tests were introduced in Guo and Chen (2016) for high-dimensional GLMs without or with nuisance regression parameters, but the p-values were obtained for testing the global hypothesis for a given set of covariates, which is different from our goal of testing the significance of individual covariates simultaneously. Portnoy (1988) studied the asymptotic behavior of the MLE for exponential families under the classical i.i.d. non-regression setting, but with diverging dimensionality. In contrast, our work under the GLM assumes the regression setting in which the design matrix X plays an important role in the asymptotic behavior of the MLE β. The validity of the asymptotic normality of the MLE was established in Portnoy (1988) under the condition of p = o(n 1/2 ), but the precise breakdown point in diverging dimensionality was not investigated therein. Another line of work is focused on generating asymptotically valid p-values when p/n converges to a fixed positive constant. For instance,  and  considered M -estimators in the linear model and showed that their variance is greater than classically predicted. Based on this result, it is possible to produce p-values by making adjustments for the inflated variance in high dimensions. Recently, Sur and Candès (2018) showed that similar adjustment is possible for the likelihood ratio test (LRT) for logistic regression. Our work differs from this line of work in two important aspects. First, our focus is on the classical p-values and their validity. Second, their results concern dimensionality that is comparable to sample size, while we aim to analyze the problem for a lower range of dimensionality and pinpoint the exact breakdown point of p-values.
The rest of the paper is organized as follows. Section 2 provides characterizations of p-values in low dimensions. We establish the nonuniformity of GLM p-values in diverging dimensions in Section 3. Section 4 presents several simulation examples verifying the theoretical phenomenon. We discuss some implications of our results in Section 5. The proofs of all the results are relegated to the Appendix.

Characterizations of P-values in Low Dimensions
To pinpoint the breakdown point of GLM p-values in diverging dimensions, we start with characterizing p-values in low dimensions. In contrast to existing work on the asymptotic distribution of the penalized MLE, our results in this section focus on the asymptotic normality of the unpenalized MLE in diverging-dimensional GLMs, which justifies the validity of conventional p-values. Although Theorems 1 and 4 to be presented in Sections 2.2 and A are in the conventional sense of relatively small p, to the best of our knowledge such results are not available in the literature before in terms of the maximum range of dimensionality p without any sparsity assumption.

Maximum likelihood estimation
For the GLM (1), the log-likelihood log f n (y; X, β) of the sample is given, up to an affine transformation, by where b(θ) = (b(θ 1 ), · · · , b(θ n )) T for θ = (θ 1 , · · · , θ n ) T ∈ R n . Denote by β = ( β 1 , · · · , β p ) T ∈ R p the MLE which is the maximizer of (2), and A well-known fact is that the n-dimensional response vector y in GLM (1) has mean vector µ(θ) and covariance matrix φΣ(θ). Clearly, the MLE β is given by the unique solution to the score equation when the design matrix X is of full column rank p. It is worth mentioning that for the linear model, the score equation (4) becomes the well-known normal equation X T y = X T Xβ which admits a closed form solution. On the other hand, equation (4) does not admit a closed form solution in general nonlinear models. This fact due to the nonlinearity of the mean function µ(·) causes the key difference between the linear and nonlinear models. In future presentations, we will occasionally use the term nonlinear GLMs to exclude the linear model from the family of GLMs when necessary.
We will present in the next two sections some sufficient conditions under which the asymptotic normality of MLE holds. In particular, Section 2.2 concerns the case of fixed design and Section A deals with the case of random design. In addition, Section 2.2 allows for general regression coefficient vector β 0 and the results extend some existing ones in the literature, while Section A assumes the global null β 0 = 0 and Gaussian random design which enable us to pinpoint the exact breakdown point of the asymptotic normality for the MLE.

Conventional p-values in low dimensions under fixed design
Recall that we condition on the design matrix X in this section. We first introduce a deviation probability bound that facilitates our technical analysis. Consider both cases of bounded responses and unbounded responses. In the latter case, assume that there exist some constants M, v 0 > 0 such that with (θ 0,1 , · · · , θ 0,n ) T = θ 0 = Xβ 0 , where β 0 = (β 0,1 , · · · , β 0,p ) T denotes the true regression coefficient vector in model (1). Then by Lv (2011, 2013), it holds that for any where ϕ(ε) = 2e −c 1 ε 2 with c 1 > 0 some constant, and ε ∈ (0, ∞) if the responses are bounded and ε ∈ (0, a 2 / a ∞ ] if the responses are unbounded. For nonlinear GLMs, the MLE β solves the nonlinear score equation (4) whose solution generally does not admit an explicit form. To address such a challenge, we construct a solution to equation (4) in an asymptotically shrinking neighborhood of β 0 that meets the MLE β thanks to the uniqueness of the solution. Specifically, define a neighborhood of β 0 as for some constant γ ∈ (0, 1/2]. Assume that p = O(n α 0 ) for some α 0 ∈ (0, γ) and let b n = o{min(n 1/2−γ √ log n, s −1 n n 2γ−α 0 −1/2 /(log n) 2 } be a diverging sequence of positive numbers, where s n is a sequence of positive numbers that will be specified in Theorem 1 below. We need some basic regularity conditions to establish the asymptotic normality of the MLE β.

Condition 1
The design matrix X satisfies with • denoting the Hadamard product and derivatives understood componentwise. Assume that max p j=1 x j ∞ < c 1/2 1 {n/(log n)} 1/2 if the responses are unbounded.

Condition 2
The eigenvalues of n −1 A n are bounded away from 0 and ∞, Conditions 1 and 2 put some basic restrictions on the design matrix X and a moment condition on the responses. For the case of linear model, bound (8) becomes (X T X) −1 ∞ = O(b n /n) and bound (9) holds automatically since b (θ) ≡ 0. Condition 2 is related to the Lyapunov condition.
Theorem 1 (Asymptotic normality) Assume that Conditions 1-2 and probability bound (6) hold. Then a) there exists a unique solution β to score equation (4) in N 0 with asymptotic probability one; b) the MLE β satisfies that for each vector u ∈ R p with u 2 = 1 and u 1 = O(s n ), and specifically for each 1 ≤ j ≤ p, where A n = X T Σ(θ 0 )X and (A −1 n ) jj denotes the jth diagonal entry of matrix A −1 n .
Theorem 1 establishes the asymptotic normality of the MLE and consequently justifies the validity of the conventional p-values in low dimensions. Note that for simplicity, we present here only the marginal asymptotic normality, and the joint asymptotic normality also holds for the projection of the MLE onto any fixed-dimensional subspace. This result can also be extended to the case of misspecified models; see, for example, Lv and Liu (2014).
As mentioned in the Introduction, the asymptotic normality was shown in Fan and Lv (2011) for nonconcave penalized MLE having intrinsic dimensionality s = o(n 1/3 ). In contrast, our result in Theorem 1 allows for the scenario of p = o(n 1/2 ) with no sparsity assumption in view of our technical conditions. In particular, we see that the conventional p-values in GLMs generally remain valid in the regime of slowly diverging dimensionality p = o(n 1/2 ).

Nonuniformity of GLM P-values in Diverging Dimensions
So far we have seen that for nonlinear GLMs, the p-values can be valid when p = o(n 1/2 ) as shown in Section 2, and can become meaningless when p ≥ n/2 as discussed in the Introduction. Apparently, there is a big gap between these two regimes of growth of dimensionality p. To provide some guidance on the practical use of p-values in nonlinear GLMs, it is of crucial importance to characterize their breakdown point. To highlight the main message with simplified technical presentation, hereafter we content ourselves with the specific case of logistic regression model for binary response. Moreover, we investigate the distributional property in (11) for the scenario of true regression coefficient vector β 0 = 0, that is, under the global null. We argue that this specific model is sufficient for our purpose because if the conventional p-values derived from MLEs fail (i.e., (11) fails) for at least one β 0 (in particular β 0 = 0), then conventional p-values are not justified. Therefore, the breakdown point for logistic regression is at least the breakdown point for general nonlinear GLMs. This argument is fundamentally different from that of proving the overall validity of conventional p-values, where one needs to prove the asymptotic normality of MLEs under general GLMs rather than any specific model.

The wild side of nonlinear regime
For the logistic regression model (1), we have b(θ) = log(1+e θ ), θ ∈ R and φ = 1. The mean vector µ(θ) and covariance matrix φΣ(θ) of the n-dimensional response vector y given by (1 + e θn ) 2 with θ = (θ 1 , · · · , θ n ) T = Xβ. In many real applications, one would like to interpret the significance of each individual covariate produced by algorithms based on the conventional asymptotic normality of the MLE as established in Theorem 1. As argued at the beginning of this section, in order to justify the validity of p-values in GLMs, the underlying theory should at least ensure that the distributional property (11) holds for logistic regression under the global null. As we will see empirically in Section 4, as the dimensionality increases, pvalues from logistic regression under the global null have a distribution that is skewed more and more toward zero. Consequently, classical hypothesis testing methods which reject the null hypothesis when p-value is less than the pre-specified level α would result in more false discoveries than the desired level of α. As a result, practitioners may simply lose the theoretical backup and the resulting decisions based on the p-values can become ineffective or even misleading. For this reason, it is important and helpful to identify the breakdown point of p-values in diverging-dimensional logistic regression model under the global null.
Characterizing the breakdown point of p-values in nonlinear GLMs is highly nontrivial and challenging. First, the nonlinearity generally causes the MLE to take no analytical form, which makes it difficult to analyze its behavior in diverging dimensions. Second, conventional probabilistic arguments for establishing the central limit theorem of MLE only enable us to see when the distributional property holds, but not exactly at what point it fails. To address these important challenges, we introduce novel geometric and probabilistic arguments presented later in the proofs of Theorems 2-3 that provide a rather delicate analysis of the MLE. In particular, our arguments unveil that the early breakdown point of p-values in nonlinear GLMs is essentially due to the nonlinearity of the mean function µ(·). This shows that p-values can behave wildly much early on in diverging dimensions when we move away from linear regression model to nonlinear regression models such as the widely applied logistic regression; see the Introduction for detailed discussions on the p-values in diverging-dimensional linear models.
Before presenting the main results, let us look at the specific case of logistic regression model under the global null. In such a scenario, it holds that θ 0 = Xβ 0 = 0 and thus Σ(θ 0 ) = 4 −1 I n , which results in In particular, we see that when n −1 X T X is close to the identity matrix I p , the asymptotic standard deviation of the jth component β j of the MLE β is close to 2n −1/2 when the asymptotic theory in (11) holds. As mentioned in the Introduction, when p ≥ n/2 the MLE can blow up with excessively large variance, a strong evidence against the distributional property in (11). In fact, one can also observe inflated variance of the MLE relative to what is predicted by the asymptotic theory in (11) even when the dimensionality p grows at a slower rate with sample size n. As a consequence, the conventional p-values given by algorithms according to property (11) can be much biased toward zero and thus produce more significant discoveries than the truth. Such a breakdown of conventional p-values is delineated clearly in the simulation examples presented in Section 4.

Main results
We now present the formal results on the invalidity of GLM p-values in diverging dimensions.
Theorem 2 (Uniform orthonormal design) 1 Assume that n −1/2 X is uniformly distributed on the Stiefel manifold V p (R n ) consisting of all n × p orthonormal matrices. Then for the logistic regression model under the global null, the asymptotic normality of the MLE established in (11) fails to hold when p ∼ n 2/3 , where ∼ stands for asymptotic order.
Theorem 3 (Correlated Gaussian design) Assume that X ∼ N (0, I n ⊗ Σ) with covariance matrix Σ nonsingular. Then for the logistic regression model under the global null, the same conclusion as in Theorem 2 holds.
Theorem 4 in Appendix A states that under the global null in GLM with Gaussian design, the p-value based on the MLE remains valid as long as the dimensionality p diverges with n at a slower rate than n 2/3 . This together with Theorems 2 and 3 shows that under the global null, the exact breakdown point for the uniformity of p-value is n 2/3 . We acknowledge that these results are mainly for theoretical interests because in practice one cannot check precisely whether the global null assumption holds or not. However, these results clearly suggest that in GLM with diverging dimensionality, one needs to be very cautious when using p-values based on the MLE.
The key ingredients of our new geometric and probabilistic arguments are demonstrated in the proof of Theorem 2 in Section B.3. The assumption that the rescaled random design matrix n −1/2 X has the Haar measure on the Stiefel manifold V p (R n ) greatly facilitates our technical analysis. The major theoretical finding is that the nonlinearity of the mean function µ(·) can be negligible in determining the asymptotic distribution of MLE as given in (11) when the dimensionality p grows at a slower rate than n 2/3 , but such nonlinearity can become dominant and deform the conventional asymptotic normality when p grows at rate n 2/3 or faster. See the last paragraph of Section B.3 for more detailed in-depth discussions on such an interesting phenomenon. Furthermore, the global null assumption is a crucial component of our geometric and probabilistic argument. The global null assumption along with the distributional assumption on the design matrix ensures the symmetry property of the MLE and the useful fact that the MLE can be asymptotically independent of the random design matrix. In the absence of such an assumption, we may suspect that p-values of the noise variables can be affected by the signal variables due to asymmetry. Indeed, our simulation study in Section 4 reveals that as the number of signal variables increases, the breakdown point of the p-values occurs even earlier.
Theorem 3 further establishes that the invalidity of GLM p-values in high dimensions beyond the scenario of orthonormal design matrices considered in Theorem 2. The breakdown of the conventional p-values occurs regardless of the correlation structure of the covariates.
Our theoretical derivations detailed in the Appendix also suggest that the conventional p-values in nonlinear GLMs can generally fail to be valid when p ∼ n α 0 with α 0 ranging between 1/2 and 2/3, which differs significantly from the phenomenon for linear models as discussed in the Introduction. The special feature of logistic regression model that the variance function b (θ) takes the maximum value 1/4 at natural parameter θ = 0 leads to a higher transition point of p ∼ n α 0 with α 0 = 2/3 for the case of global null β 0 = 0.

Numerical Studies
We now investigate the breakdown point of p-values for nonlinear GLMs in diverging dimensions as predicted by our major theoretical results in Section 3 with several simulation examples. Indeed, these theoretical results are well supported by the numerical studies.

Simulation examples
Following Theorems 2-3 in Section 3, we consider three examples of the logistic regression model (1). The response vector y = (y 1 , · · · , y n ) T has independent components and each y i has Bernoulli distribution with parameter e θ i /(1 + e θ i ), where θ = (θ 1 , · · · , θ n ) T = Xβ 0 . In example 1, we generate the n × p design matrix X = (x 1 , · · · , x p ) such that n −1/2 X is uniformly distributed on the Stiefel manifold V p (R n ) as in Theorem 2, while examples 2 and 3 assume that X ∼ N (0, I n ⊗ Σ) with covariance matrix Σ as in Theorem 3. In particular, we choose Σ = (ρ |j−k| ) 1≤j,k≤p with ρ = 0, 0.5, and 0.8 to reflect low, moderate, and high correlation levels among the covariates. Moreover, examples 1 and 2 assume the global null model with β 0 = 0 following our theoretical results, whereas example 3 allows sparsity s = β 0 0 to vary.
To examine the asymptotic results we set the sample size n = 1000. In each example, we consider a spectrum of dimensionality p with varying rate of growth with sample size n. As mentioned in the Introduction, the phase transition of perfect hyperplane separation happens at the point p/n = 1/2. Recall that Theorems 2-3 establish that the conventional GLM p-values can become invalid when p ∼ n 2/3 . We set p = [n α 0 ] with α 0 in the grid {2/3 − 4δ, · · · , 2/3 − δ, 2/3, 2/3 + δ, · · · , 2/3 + 4δ, (log(n) − log(2))/ log(n)} for δ = 0.05. For example 3, we pick s signals uniformly at random among all but the first components, where a random half of them are chosen as 3 and the other half are set as −3.
The goal of the simulation examples is to investigate empirically when the conventional GLM p-values could break down in diverging dimensions. When the asymptotic theory for the MLE in (11) holds, the conventional p-values would be valid and distributed uniformly on the interval [0, 1] under the null hypothesis. Note that the first covariate x 1 is a null variable in each simulation example. Thus in each replication, we calculate the conventional p-value for testing the null hypothesis H 0 : β 0,1 = 0. To check the validity of these p-values, we further test their uniformity.
For each simulation example, we first calculate the p-values for a total of 1, 000 replications as described above and then test the uniformity of these 1, 000 p-values using, for example, the Kolmogorov-Smirnov (KS) test (Kolmogorov, 1933;Smirnov, 1948) and the Anderson-Darling (AD) test Darling, 1952, 1954). We repeat this procedure 100 times to obtain a final set of 100 new p-values from each of these two uniformity tests. Specifically, the KS and AD test statistics for testing the uniformity on [0, 1] are defined as

Testing results
For each simulation example, we apply both KS and AD tests to verify the asymptotic theory for the MLE in (11) by testing the uniformity of conventional p-values at significance level 0.05. As mentioned in Section 4.1, we end up with two sets of 100 new p-values from the KS and AD tests. Figures 1-3  regression model under global null that the conventional p-values break down when p ∼ n α 0 with α 0 = 2/3. Figure 3 for example 3 examines the breakdown point of p-values with varying sparsity s. It is interesting to see that the breakdown point shifts even earlier when s increases as suggested in the discussions in Section 3.2. The results from the AD test are similar so we present only the results from the KS test for simplicity.
To gain further insights into the nonuniformity of the null p-values, we next provide an additional figure in the setting of simulation example 1. Specifically, in Figure 4 we present the histograms of the 1,000 null p-values from the first simulation repetition (out of 100) for each value of α 0 . It is seen that as the dimensionality increases (i.e., α 0 increases), the null p-values have a distribution that is skewed more and more toward zero, which is prone to produce more false discoveries if these p-values are used naively in classical hypothesis testing methods.
To further demonstrate the severity of the problem, we estimate the probability of making type I error at significance level a, as the fraction of p-values below a. The means q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  Table 1 for a = 0.05 and 0.1. When the null p-values are distributed uniformly, the probabilities of making type I error should all be close to the target level a. However, Table 1 shows that when the growth rate of dimensionality α 0 approaches or exceeds 2/3, these probabilities can be much larger than a, which again supports our theoretical findings. Also it is seen that when α 0 is close to but still smaller than 2/3, the averages of estimated probabilities exceed slightly a, which could be the effect of finite sample size.

Discussions
In this paper we have provided characterizations of p-values in nonlinear GLMs with diverging dimensionality. The major findings are that the conventional p-values can remain valid when p = o(n 1/2 ), but can become invalid much earlier in nonlinear models than in linear models, where the latter case can allow for p = o(n). In particular, our theoretical results pinpoint the breakdown point of p ∼ n 2/3 for p-values in diverging-dimensional logistic regression model under global null with uniform orthonormal design and correlated Gaussian design, as evidenced in the numerical results. It would be interesting to investigate such a phenomenon for more general class of random design matrices. The problem of identifying the breakdown point of p-values becomes even more complicated and challenging when we move away from the setting of global null. Our technical analysis suggests that the breakdown point p ∼ n α 0 can shift even earlier with α 0 ranging between 1/2 and 2/3. But the exact breakdown point can depend upon the number of signals s, the signal magnitude, and the correlation structure among the covariates in a rather complicated fashion. Thus more delicate mathematical analysis is needed to obtain the exact relationship. We leave such a problem for future investigation. Moving beyond the GLM setting will further complicate the theoretical analysis.
As we routinely produce p-values using algorithms, the phenomenon of nonuniformity of p-values occurring early in diverging dimensions unveiled in the paper poses useful cautions to researchers and practitioners when making decisions in real applications using results from p-value based methods. For instance, when testing the joint significance of covariates in diverging-dimensional nonlinear models, the effective sample size requirement should be checked before interpreting the testing results. Indeed, statistical inference in general highdimensional nonlinear models is particularly challenging since obtaining accurate p-values is generally not easy. One possible route is to bypass the use of p-values in certain tasks including the false discovery rate (FDR) control; see, for example, Barber and Candès (2015); ;  for some initial efforts made along this line.

Acknowledgments
This work was supported by NIH Grant 1R01GM131407-01, NSF CAREER Award DMS-1150318, a grant from the Simons Foundation, and Adobe Data Science Research Award. The first and last authors sincerely thank Emmanuel Candès for helpful discussions on this topic. The authors would like to thank the Associate Editor and referees for their valuable comments that helped improve the article substantially.

Appendix A. Conventional P-values in Low Dimensions under Random Design
Under the specific assumption of Gaussian design and global null β 0 = 0, we can show that the asymptotic normality of MLE continues to hold without previous Conditions 1-2.
Theorem 4 shows that the conclusions of Theorem 1 continue to hold for the case of random design and global null with the major difference that the dimensionality can be pushed as far as p ∼ n 2/3 . The main reasons for presenting Theorem 4 under Gausssian design are twofold. First, Gaussian design is a widely used assumption in the literature. Second, our results on the nonuniformity of GLM p-values in diverging dimensions use geometric and probabilistic arguments which require random design setting; see Section 3 for more details. To contrast more accurately the two regimes and maintain self-contained theory, we have chosen to present Theorem 4 under Gaussian design. On the other hand, we would like to point out that Theorem 4 is not for practitioners who want to justify the usage of classical p-values. The global null assumption of β 0 = 0 restricts the validity of Theorem 4 in many practical scenarios.

Appendix B. Proofs of Main Results
We provide the detailed proofs of Theorems 1-3 in this Appendix.

B.1. Proof of Theorem 1
To ease the presentation, we split the proof into two parts, where the first part locates the MLE β in an asymptotically shrinking neighborhood N 0 of the true regression coefficient vector β 0 with significant probability and the second part further establishes its asymptotic normality.
Part 2: Conventional asymptotic normality of the MLE β. Fix any 1 ≤ j ≤ p. In light of (16), we have β − β 0 = A −1 n (ξ − r), which results in with e j ∈ R p having one for the jth component and zero otherwise. Note that since the smallest and largest eigenvalues of n −1 A n are bounded away from 0 and ∞ by Condition 2, it is easy to show that (A −1 n ) −1/2 jj is of exact order n 1/2 . In view of (17), it holds on the event E defined in (12) that since b n = o{n 2γ−α 0 −1/2 /(log n) 2 } by assumption. This leads to It remains to consider the term (A −1 n ) Clearly, the n random variables η i 's are independent with mean 0 and It follows from Condition 2 and the Cauchy-Schwarz inequality that Thus an application of Lyapunov's theorem yields By Slutsky's lemma, we see from (20)-(22) that showing the asymptotic normality of each component β j of the MLE β.
We further establish the asymptotic normality for the one-dimensional projections of the MLE β. Fix an arbitrary vector u ∈ R p with u 2 = 1 satisfying the L 1 sparsity bound u 1 = O(s n ). In light of (16), we have β − β 0 = A −1 n (ξ − r), which results in Note that since the smallest and largest eigenvalues of n −1 A n are bounded away from 0 and ∞ by Condition 2, it is easy to show that (u T A −1 n u) −1/2 is of exact order n 1/2 . In view of (17), it holds on the event E defined in (12) that since b n = o{s −1 n n 2γ−α 0 −1/2 /(log n) 2 } by assumption. This leads to since u 1 = O(s n ) by assumption. It remains to consider the term (u T A −1 Clearly, the n random variables η i 's are independent with mean 0 and It follows from Condition 2 and the Cauchy-Schwarz inequality that Thus an application of Lyapunov's theorem yields By Slutsky's lemma, we see from (23)-(25) that showing the asymptotic normality of any L 1 -sparse one-dimensional projection u T β of the MLE β. This completes the proof of Theorem 1.
Thus, the MLE must fall into the region N 0 following the similar arguments in Theorem 1. Next, we show the componentwise asymptotic normality of the MLE β. By equation where . By Lemma 13 and Equation (32), both n 1/2 (b (0)) −1 T and n 1/2 e T j [b (0)X T X] −1 (r+s+t) converges to zero in probability. So, it is enough to consider the first summand in (33). Now, we show that n −1/2 e T j X T [y − µ 0 ] is asymptotically normal. In fact, we can write e T j X T [y−µ 0 ] = n i=1 x ij y i where each summand x ij y i is independent over i and has variance φb (0). Moreover, n i=1 E|x ij y i | 3 = O(n) since |x ij | 3 and |y i | 3 are independent and finite mean. So, we apply Lyapunov's theorem to ). Finally, we know that b (0)n(A −1 n ) jj → 1 in probability from the remark in Theorem 1. Thus, Slutsky's lemma yields This completes the proof of the theorem.
Lemma 5 Assume that the components of y−µ 0 are uniform sub-Gaussians. That is, there exist a positive constant C such that P (|(y − µ 0 ) i | > t) ≤ C exp{−Ct 2 } for all 1 ≤ i ≤ n. Then, it holds that, for some positive constant c 2 , with asymptotic probability 1 − o(p −a ).
Proof We prove the result by conditioning on X. Let E = n −1 X T X − I p . Then by matrix inversion, Thus, it follows that In the rest of the proof, we will bound η 1 , η 2 and η 3 . Part 1: Bound of η 1 . First, it is easy to see that We observe that each summand x ij (y − µ 0 ) i is the product of two subgaussian random variables, and so satisfies P (|x ij (y − µ 0 ) i | > t) ≤ C exp(−Ct) for some constant C > 0 by Lemma 1 in Fan et al. (2016). Moreover, E[x ij (y − µ 0 ) i ] = 0 since x ij and (y − µ 0 ) i are independent and have zero mean. Thus, we can use Lemma 9 by setting W ij = x ij (y − µ 0 ) i and α = 1. So, we get with probability 1 − O(p −c ) for some positive constants c and c 2 . Part 2: Bound of η 2 . Now, we study η 2 = n −2 X T XX T (y − µ 0 ) ∞ . Let z k be the k-th column of X, that is z k = Xe k . Direct calculations yield By Lemma 14, we have max k z k 2 = max k z k 2 2 ≤ O p ( √ n). Therefore, by using the fact that y Combining (36)-(38) yields Part 3: Bound of η 3 . Finally, we study η 3 . We observe that (35) shows that n −1 X T (y − µ 0 ) ∞ = O( n −1 log p) with probability 1 − O(p −c ). Putting these facts together, we obtain where we use p = O(n α 0 ) with α 0 ∈ [0, 2/3). Combining equations (35), (39), and (40), we obtain that with probability at least 1 − o(p −a ), (X T X) −1 X T (y − µ 0 ) ∞ ≤ c n −1 log n.
Lemma 6 Under the assumptions of Theorem 4, Proof Let E = n −1 X T X − I p . Then, E 2 ≤ C(p/n) 1/2 for some constant C with probability 1 − O(p −c ) by Theorem 4.6.1 in Vershynin (2016). Furthermore, by matrix inversion, we get Now, we take the norm and use triangle inequalities to get where we use the fact that p/n is bounded by a constant less than 1.
Lemma 7 In the same setting as Lemma 6, if E = n −1 X T X−I p , then Proof Again, we use that E 2 ≤ C(p/n) 1/2 for some constant C with probability 1 − O(p −c ). By similar calculations as in Lemma 6, we deduce Lemma 8 Let W j be nonnegative random variables for 1 ≤ j ≤ p that are not necessarily independent. If P (W j > t) ≤ C 1 exp(−C 2 a n t 2 ) for some constants C 1 and C 2 and for some sequence a n , then for any c > 0, max 1≤j≤p W j ≤ ((c + 1)/C 2 ) 1/2 a −1/2 n (log p) 1/2 with probability at least 1 − O(p −c ).
Proof Using union bound, we get Taking t = a −1/2 n (log p) 1/2 ((c + 1)/C 2 ) 1/2 concludes the proof since then Lemma 9 Let W ij be random variables which are independent over the index i. Assume that there are constants C 1 and C 2 such that P (|W ij | > t) ≤ C 1 exp(−C 2 t α ) with 0 < α ≤ 1.
for some positive constants c and C. Fan et al. (2016) where C 3 and C 4 are some positive constants which only depend on C 1 and C 2 . This probability bound shows that the assumption of Lemma 8 holds with a n = n α . Thus, Lemma 8 finishes the proof.
Proof First, observe that for some constant C, Moreover, the summands x ij x T i β β 2 3 are independent over i and they satisfy the probability bound P (|x ij x T i β β 2 3 | > t) ≤ C exp(−Ct 1/2 ) by Lemma 1 of Fan et al. (2016). Thus, by Lemma 9, we obtain = O(n −1/4 (log p) 1/2 ). Now, we calculate the expected value of the summand x ij x T i β β 2

3
. We decompose x T i β as x ij β j + x T i,−j β −j where x i,−j and β −j are the vectors x i and β whose jth entry is removed. We use the independence of x i,−j and x ij and get Finally, we can combine the result of Lemma 9 and the expected value of x ij x T i β β 2
Under the assumptions of Theorem 4, we have Proof Since X and y are independent, expectation of T is clearly zero. Then, we consider the variance of T . To this end, we condition on X. We can calculate the conditional variance of T as follows Now, we can obtain the unconditional variance using the law of total variance.

var[T ] = E[var[T |X]] + var[E[T |X]]
Thus, using Lemma 7, we can show that var[T ] = o(n −1 ). Finally, we use Chebyshev's inequality P (|T | > n −1/2 ) ≤ nvar[T ] = o(1). So, we conclude that T = o p (n −1/2 ) Lemma 14 Let x ij be standard normal random variables for 1 ≤ i ≤ n and 1 ≤ j ≤ p. Then, max 1≤j≤p n i=1 x 2 ij ≤ n + O(n 1/2 (log p) 1/2 ) with probability 1 − O(p −c ) for some positive constant c. Consequently, when log p = O(n α ) for some 0 < α ≤ 1, we have max 1≤j≤p n i=1 x 2 ij = O(n), for large enough n with probability 1 − O(p −c ). Proof Since x ij is a standard normal variable, x 2 ij is subexponential random variable whose mean is 1. So, Lemma 9 entails that The remaining analysis focuses on the score equation (4) which is solved exactly by the MLE β, that is, which leads to Let us first consider the random variable ξ defined in (51). Note that 2[y − µ(0)] has independent and identically distributed (i.i.d.) components each taking value 1 or −1 with equal probability 1/2, and is independent of X. Thus since n −1/2 X is uniformly distributed on the Stiefel manifold V p (R n ), it is easy to see that where 1 ∈ R n is a vector with all components being one. Using similar arguments as before, we can show that ξ has a spherical distribution on R p . Thus the joint distribution of ξ is determined completely by the marginal distribution of ξ. For each 1 ≤ j ≤ p, denote by ξ j the jth component of ξ = 2 −1 n −1/2 X T 1 using the distributional representation in (52). Let X = (x 1 , · · · , x p ) with each x j ∈ R n . Then we have ξ j = 2 −1 n −1/2 x T j 1 d = 2 −1 (n 1/2 / x j 2 )n −1/2 x T j 1, where x j ∼ N (0, 4 −1 I n ). It follows from (53) and the concentration phenomenon of Gaussian measures that each ξ j is asymptotically close to N (0, 4 −1 ) and thus consequently ξ is asymptotically close to N (0, 4 −1 I p ). A key fact (i) for the finite-sample distribution of ξ is that the standard deviation of each component ξ j converges to 1/2 at rate O P (n −1/2 ) that does not depend upon the dimensionality p at all. We now turn our attention to the second term η defined in (51). In view of (49) and the fact that n −1/2 X is uniformly distributed on the Stiefel manifold V p (R n ), we can show that with significant probability, X β ∞ ≤ o(1) for p ∼ n α 0 with α 0 < 1. The uniform bound in (54) enables us to apply the mean value theorem for the vector-valued function η around β 0 = 0, which results in = 4 −1 n 1/2 β + r since n −1/2 X is assumed to be orthonormal, where r = n −1/2 X T 1 0 Σ(tX β) − 4 −1 I n dt X β.
Here, the remainder term r = (r 1 , · · · , r p ) T ∈ R p is stochastic and each component r j is generally of order O P {p 1/2 n −1/2 } in light of (49) when the true model may deviate from the global null case of β 0 = 0.
Since our focus in this theorem is the logistic regression model under the global null, we can in fact claim that each component r j is generally of order O P {pn −1 }, which is a better rate of convergence than the one mentioned above thanks to the assumption of β 0 = 0. To prove this claim, note that the variance function b (θ) is symmetric in θ ∈ R and takes the maximum value 1/4 at θ = 0. Thus in view of (54), we can show that with significant probability, 4 −1 I n − Σ(tX β) ≥ cdiag{(tX β) • (tX β)} = ct 2 diag{(X β) • (X β)} for all t ∈ [0, 1], where c > 0 is some constant and ≥ stands for the inequality for positive semidefinite matrices. Moreover, it follows from (49) and the fact that n −1/2 X is uniformly distributed on the Stiefel manifold V p (R n ) that with significant probability, all the n components of X β are concentrated in the order of p 1/2 n −1/2 . This result along with (57) and the fact that n −1 X T X = I p entails that with significant probability, ≥ n −1/2 X T 1 0 c * t 2 pn −1 dt X = 3 −1 c * pn −3/2 X T X = 3 −1 c * pn −1/2 I p , where c * > 0 is some constant. Thus combining (56), (58), and (49) proves the above claim. We make two important observations about the remainder term r in (55). First, r has a spherical distribution on R p . This is because by (55) and (51) it holds that r = η − 4 −1 n 1/2 β = ξ − 4 −1 n 1/2 β, which has a spherical distribution on R p . Thus the joint distribution of r is determined completely by the marginal distribution of r. Second, for the nonlinear setting of logistic regression model, the appearance of the remainder term r in (55) is due solely to the nonlinearity of the mean function µ(·), and we have shown that each component r j can indeed achieve the worst-case order pn −1 in probability. For each 1 ≤ j ≤ p, denote by η j the jth component of η. Then in view of (49) and (55), a key fact (ii) for the finite-sample distribution of η is that the standard deviation of each component η j converges to 1/2 at rate O P {pn −1 } that generally does depend upon the dimensionality p.
Finally, we are ready to compare the two random variables ξ and η on the two sides of equation (51). Since equation (51) is a distributional identity in R p , naturally the square root of the sum of varξ j 's and the square root of the sum of varη j 's are expected to converge to the common value 2 −1 p 1/2 at rates that are asymptotically negligible. However, the former has rate p 1/2 O P (n −1/2 ) = O P {p 1/2 n −1/2 }, whereas the latter has rate p 1/2 O P {pn −1 } = O P {p 3/2 n −1 }. A key consequence is that when p ∼ n α 0 for some constant 2/3 ≤ α 0 < 1, there is a profound difference between the two asymptotic rates in that the former rate is O P {n −(1−α 0 )/2 } = o P (1), while the latter rate becomes O P {n 3α 0 /2−1 } which is now asymptotically diverging or nonvanishing. Such an intrinsic asymptotic difference is, however, prohibited by the distributional identity (51) in R p , which results in a contradiction. Therefore, we have now argued that assumption (A) we started with for probability one, n −1/2 X T 1 0 4 −1 I n − Σ(tX β) dt X (63) ≥ n −1/2 X T 1 0 c * t 2 pn −1 dt X = 3 −1 c * pn −3/2 X T X → 3 −1 c * pn −1/2 I p , where c * > 0 is some constant. This completes the proof of Theorem 3.